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The association of ions in electrolyte solutions at very low concentration and low temperature is 
studied using computer simulations and quasi-chemical ion-pairing theory. The specific case of the 
restricted primitive model (charged hard spheres) is considered. Specialised simulation techniques 
are employed that lead to efficient sampling of the arrangements and distributions of clusters and 
free ions, even at conditions corresponding to nanomolar solutions of simple salts in solvents with 
dielectric constants in the range 5-10, as used in recent experimental work on charged-coUoid sus- 
pensions. A direct comparison is effected between theory and simulation using a variety of clustering 
criteria and theoretical approximations. It is shown that conventional distance-based cluster criteria 
can give erroneous results. A reliable set of theoretical and simulation estimators for the degree of 
association is proposed. The ion-pairing theory is then compared to experimental results for salt 
solutions in low-polarity solvents. The agreement is excellent, and on this basis some calculations 
are made for the screening lengths which will figure in the treatment of colloid-colloid interactions 
in such solutions. The accord with available experimental results is complete. 



I. INTRODUCTION 

The control of chargcd-colloid suspensions with added 
salt is a linchpin of soft condensed matter science. The 
physical principles that govern the net interactions be- 
tween like-charged colloids in aqueous electrolyte solu- 
tions were laid down more than half a century agoii 
Central to charge stabilization are the formation of the 
electrical double layer and the phenomenon of screen- 
ing over distances comparable to the Dcbye length. The 
classic Derjaguin-Landau-Verwey-Overbeek and Debye- 
Hiickel (DH) theories apply well to high-polarity solvents, 
such as water, where the electrostatic interactions be- 
tween salt ions and counterions are strongly screened di- 
electrically. In particular, under normal conditions and 
with simple salts, cation-anion pairing is rather insignifi- 
cant and complete ion dissociation can be assumed when 
describing the screening effect. 

Charged colloidal suspensions in low dielectric constant 
solvents are now of experimental interest i^"— Such sys- 
tems involve whole new regimes of low ion concentrations 
~ nM, and strong electrostatic interactions between ions 
(added salt and counterions). The situation with regard 
to screening now changes because the reduction in di- 
electric screening compared to that in water can only en- 
hance cation-anion association and promote the forma- 
tion of a significant number of so-called Bjerrum pairs. 
This leads to a reduction in free ions and a concomitant 
increase in the screening length. The effects of ion pair- 
ing on the screening of colloidal interactions have been 
cxplorcdJ^"— 

Of course, ion pairing is not a new phenomenon, and its 
effects on the thermodynamics and dynamical properties 
can be highly pronouncedJ^ Already in 1926, Bjerrum 
described his eponymous pairs within a quasi-chemical 
ion-pairing equilibrium, and suggested removing them 
from the effective free-ion concentration when performing 



DH-like calculations on electrolyte solutions Bjerrum's 
theory has been thoroughly tested against experimental 
data for solutions with moderately low salt concentra- 
tions of Cs > 10~^ M in solvents with dielectric constants 
in the range 2 < e < 80J^ One of the most dramatic 
manifestations of ion pairing is in the phase separation 
of ionic fluidsfiiii^ where the low-concentration 'vapour' 
phase has such a high degree of ion association that the 
conventional DH theory has to be extended to include 
ion-dipole and dipole-dipole interactions in order to give 
a good account of accurate simulation data for the coex- 
istence cnvelope i^^'^° 

From a computational perspective, the new experi- 
mental regimes of very low concentration and strong 
electrostatic interactions present some serious challenges. 
Molecular dynamics simulations of salts at nanomo- 
lar concentrations have already fallen foul of sampling 
problems J ^^'^'^ Recently, the current authors put forward 
a protocol for performing Monte Carlo (MC) simulations 
in the canonical ensemble, with novel particle moves that 
allow efficient equilibration at the extreme conditions re- 
ferred to abovej^i This opens up the opportunity to ex- 
plore the true degree of association in very low concentra- 
tion electrolyte solutions made up with low-polarity sol- 
vents. For the purposes of this exploratory study, atten- 
tion is focused on the restricted primitive model (RPM) 
of ionic fluids. The RPM is an electroneutral mixture of 
N/2 positively and iV/2 negatively charged hard spheres 
of equal diameter a and charges ibg immersed in a di- 
electric continuum with dielectric constant e and volume 
V at temperature T, with an overall ion concentration 
p = N/V. The interaction pair potential between ions i 
and j is 

{oo Tij < a 
^r.,,>a (1) 

where is the pair separation, g,; is the charge on ion 
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i, and D = A^ne^e where eo is the dielectric permittiv- 
ity of vacuum. The overall ion concentration and tem- 
perature are given in reduced units by p* = pa^ and 
T* = kBTDa/q^, respectively. The Bjerrum length is 
the distance at which the attractive cation-anion poten- 
tial is equal to —k^T, and is given by Ab = <j/T* . The 
phase behaviour of the RPM is now well known j^^i^'^ the 
vapour-liquid critical parameters are T* ~ 0.05 and ~ 
O.OSi ^^i^^ Simulations confirm that the degree of ion asso- 
ciation in the vapour phase just below p* and T* ~ 0.05 
is significant)^"— Indeed, the coexistence properties of 
fused cation-anion pairs (charged hard dumbbells) are al- 
most identical to those of the RPM,^"i2i 

In this work, the degree of ion association in the RPM 
at very low concentrations and low (near-critical) temper- 
atures is investigated. Calculations are performed down 
to a reduced ion concentration of p* = 10~^° and a re- 
duced temperature of T* = 0.04; for a monovalent salt 
with ionic diameter cr = 4 A at room temperature, these 
values correspond to a salt concentration Cg — 1 nM and a 
solvent dielectric constant e ~ 5.6. Using specialised MC 
simulations, results are obtained with which to test the 
quasi-chemical ion-pairing theory as proposed by Bjer- 
rum. This involves using a novel simulation protocol re- 
cently proposed by us;^ and a variety of methods for 
determining the degree of association. The ion-pairing 
theory is then tested against experimental data for the 
degree of ion association; to this end, recent work by Le- 
unissen and co-workers has yielded results for salts at con- 
centrations of around 10~^ M in solvents with dielectric 
constants as low as about 5.— To the best of our knowl- 
edge, this is the first time that a quantitative comparison 
has been made between theory, simulation, and experi- 
ment at such extreme conditions. On the basis of this 
comparison, the effects of ion association on the screen- 
ing of charged- colloid interactions under such conditions 
can be evaluated with some confidence. 

This article is arranged as follows. The ion-pairing the- 
ory is presented in Section [III and the simulation details 
are summarised in Section IIIII The RPM simulation re- 
sults are given in Section IIV A[ and an analysis of ex- 
perimental data is presented in Section HVBI Section IVl 
concludes the paper. 



II. THEORY 

To describe ion association, consider the quasi-chemical 
equilibrium 



The choice of the cutoff radius is to be discussed below. 
At equilibrium p,± = /t+ + which leads to 



cation-anion pair ^ cation + anion. 



(2) 



In terms of the degree of association a, the concentra- 
tion of cation-anion pairs is p± = Oip/2 and those of 
the cations and anions are p+ = p^ = (1 — a)p/2. 
Considering the mixture of cations, anions, and cation- 
anion pairs to be ideal, the chemical potentials are p,± = 
kBT\n{apK\k^_ /2K), p+ = fceTln [(1 - a)pA3 /2], and 
/i_ = /cbTIu [(1 — a)/9A'l/2], where A+ and A_ are the 
de Broglie thermal wavelengths of the cations and anions, 
respectively. Here K is the configurational integral of a 
pair, which plays the role of an equilibrium constant: 



Kp 
~2~' 



Solving for a yields 
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a = l- — {^T+2Kp-l 
Kp \ 



(4) 



(5) 



Considering the phase diagram of the RPM in the 
concentration-temperature plane, a sensible dividing line 
between 'associated' and 'dissociated' regimes is the locus 
of points defined by a = ^ , or alternatively 

Kp = 4. (6) 

All that remains now is to determine the equilibrium 
constant K. The primary problem is that the integral 
in Eq. ^ does not converge for Tc — J- oo, and so in 
the conventional treatment, an appropriate finite upper 
limit for the integral needs to be identified. One choice 
for Tc is the Bjerrum length Ab, on the basis that the 
separation between ions in a pair should be such that 
the interaction energy is greater in magnitude than k^T. 
An alternative, and more conventional, choice is to set 
fc = Ab/2 corresponding to the minimum of the inte- 
grand r'^exp{XB/r) in Eq. ([3]). It has long been recog- 
nised, however, that the precise choice of Tc is unimpor- 
tant (see section 925 of Ref. IsgI ). at least at low temper- 
atures; this will be emphasised in the results of the cur- 
rent work. An approximate closed-form expression for 
K valid at low T* can be obtained by noting that in a 
cation-anion pair, the separation r should not be much 
more than cr. Writing r = a + 6r leads to the limiting 
behaviour a/r sa 1 — Sr/a = 2 — r/a. Substituting this 
in to Eq. ([3]) and performing the integral with rc = oo 
yields 



[T* + 2{T*f + 2{T*f 



(7) 



K = 4:TT 



cxp ^— ^^dr. 



(3) 



Equation (O is possibly the most simple low-temperature 
result, and was inspired by a similar approximation for 
the two-particle partition function of dipolar hard spheres 
presented by Jordan;^ the range of validity is limited by 
an unphysical minimum in K at T* ~ 0.54. Levin and 
Fisher have summarised several more accurate closed- 
form expressions)^ Finally, it is acknowledged that Ebel- 
ing's alternative expression for which reproduces the 
correct equation of state for the RPM up to terms of order 
and is therefore a more rigorous choice,— gives essen- 
tially identical results to the Bjerrum-length prescriptions 
employed here. 



III. SIMULATION METHODS 

Conventional MC simulations of associating fiuids at 
very low concentrations can fail due to insufficient sam- 
pling of the most significant arrangements and spatial dis- 
tributions of clusters.— On the one hand, during a typical 
length run using single-particle moves, isolated particles 
in very dilute systems may never come within suflacient 
proximity of other particles to associate. On the other 
hand, particles already within clusters may not be able 
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to detach due to it being a rare event. In an effort to 
eiiminate tfiese problems, the authors recently proposed 
an efficient MC protocol for simulating the RPM at very 
low concentrations and low temperatures where ion as- 
sociation is expected to be significant<2i The simulations 
are conducted within the canonical (NVT) ensemble us- 
ing a cubic box of side L = V^^^ with periodic bound- 
ary conditions. The long-range coulombic interactions 
are handled using the Ewald summation with conducting 
boundary conditions4S Various types of MC moves are 
attempted: normal single-particle moves with displace- 
ments chosen randomly from either a narrow interval 
(with a width adjusted to give an acceptance rate of 40%) 
or a broad interval (spanning the range —L/2 to L/2)] 
cluster moves with displacements chosen randomly from 
narrow and broad intervals, as before; and 'cluster forma- 
tion/breakage' (CFB) moves, each of which involves mov- 
ing a second ion within a sphere of radius A centered on a 
randomly chosen first ion. This last move offers possibil- 
ities for bringing together two randomly selected isolated 
ions in to association, and for prising two clustered ions 
apart. Full details of the simulation protocol are reported 
in Ref.[2l|. The main control parameters are the radius A, 
and the various proportions of single-particle and cluster 
moves, small and large displacements, and CFB moves. 
On the basis of earlier work;^ the present simulations are 
performed with 70% small single-particle displacements, 
10% large single-particle displacements, 5% small cluster 
displacements, 5% large cluster displacements, and 10% 
CFB moves. The CFB radius was set to A = L/4: in all 
cases. These parameters were shown in Ref. [2l| to give 
rapid convergence to the apparent equilibrium state. In 
all cases, the system is made up of = 256 ions, and run 
lengths consist of 10^-10^ MC moves per ion, depending 
on density and temperature. 

Conventionally, clusters in fluids are identifled using 
some kind of pairwise distance^i or energy-based crite- 
rion; the latter are useful for anisotropic potentials, where 
not only the distance but also the orientation have to be 
favorable for association to occur. In the present case, 
a distance-based criterion suffices; two particles are con- 
sidered associated if their separation is less than some 
cutoff distance Vc- In their comprehensive study of ion 
association in the vapour phase of the RPM near coex- 
istence, Caillol and Weis showed that the cluster distri- 
bution is basically independent of criteria in the range 
1.8(7 < Tc < 2.2fT4i Allahyarov et al. use = Ab in 
Ref. m and Tc = 3a in Ref. [l^ In this work, three 
different distance criteria were employed: rc ~ Ab and 
fc = Ab/2 are obvious candidates, for the reasons out- 
lined in Section |lll and rc = 2a, in line with earlier 
studies4i Using these criteria, a is the proportion of ions 
clustered with at least one other ion. 

In addition, outlined here is a method of estimating the 
degree of association a from simulation data without hav- 
ing to specify a cluster criterion. Consider the nearest- 
neighbour cation-anion distribution function p(r), reflect- 
ing the distance between an ion and its nearest neighbour 
of opposite charge. If a cation is dissociated, then the 
nearest-neighbour anion is remote (due to the low den- 
sities of interest here) and to a first approximation can 
be assumed completely uncorrelated with the cation. The 
probability of an anion being at a distance between r and 
r + dr from the cation, and the remaining N/2 — 1 anions 



being at least as far away, is 

Pd(r)dr = -X— (^1-^j 

« 27r/cir^ exp (-|7r/9r^)dr (8) 

where the subscript 'd' denotes 'dissociated'. Note that 
Pd{f) is normalised and shows a peak at rg = (I/tt/?)^/'^. 
If the nearest-neighbour anion is associated with the 

cation, then the radial distribution function (r) will 

be peaked near r = a, signalling very strong, short-range 
correlations which are not amenable to an accurate theo- 
retical treatment; the corresponding function for associ- 
ated ('a') cations, Pa{r), is not easy to predict. For the 
RPM the arguments above apply in exactly the same way 
to anions. If the proportion of associated ions is a, and 
that of dissociated ions is (1 — a), then the total p{r) will 
be given by 

p{r) = api,{r) + {I - a)pd{r). (9) 

This function can be obtained directly from simulations 
and, in principle, fitting Eq. ([9]) to simulation results 
yields the degree of association without having to specify 
a cluster criterion. In practice, and without a reliable ex- 
pression for Pii{r), (1 — a)p(i{r) is fitted to p{r) over the 
range r > rg, where apa(^) makes no significant contri- 
bution: 

p{r) ~ (1 - a) [2npr'^ exp (-fTrpr^)] r > vq. (10) 

IV. RESULTS 
A. Restricted primitive model 

Figure [T] shows the degree of association a as a func- 
tion of reduced ion density p* along several isotherms, 
T* = 0.04, 0.05, 0.06, and 0.07. RecaU that the critical 
temperature of the RPM is T* ~ 0.05. Four sets of simu- 
lation data are shown, corresponding to different cluster- 
ing criteria. The simulation data in Figs. [Ila)-(c) were 
obtained using distance-based criteria of Tc = Ab, Ab/2, 
and 2a; the data in Fig. [TJd) were obtained from fits to 
p{r). Also included in the figures are the theoretical pre- 
dictions of Eq. ([5]) with K evaluated numerically using 
Eq. ^ and rc = Ab/2, and with the asymptotic expres- 
sion in Eq. ([7]). On the scale of these plots, curves with 
fc = Ab are indistinguishable from those with Tc = Ab/2 
and so they are omitted. 

The first impression given by Fig. [T] is that there is 
very good overall agreement between the simulation re- 
sults and the theoretical predictions. A close inspection 
of Figs. [Ha)-(c) shows that, in simulations, the distance 
criteria rc = Ab and Ab/2 give slightly poorer results for 
the degree of association; looking at Figs. [TJa) and (b), 
the data along the higher temperature isotherms vary 
a little too sharply and saturate at a = 1 prematurely 
as the density is increased. Figure [TJc) shows that the 
fixed-distance cut-off of Tc = 2a provides a more realistic 
variation with density, refiecting a strong association of 
ions in clusters close to contact. 

Further insights are afforded by simulation measure- 
ments of the nearest-neighbour cation-anion distribution 
function, p{r). Two examples from the T* = 0.05 
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(a) r = X, 



(b) r = X / 2 
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FIG. 1: The degree of association a along isotherms. The 
curves are the theoretical predictions of Eq. ^ using K com- 
puted with a cut-off rc = Ab/2 (solid lines), and using the 
asymptotic result in Eq. (O for K (dashed lines). The points 
are the simulation results computed using various criteria: (a) 
distance criterion with Vc = Ab; (b) distance criterion with 
Vc — Ab/2; (c) distance criterion with rc — 2a; (d) fitting 
Eq. dlOj to p{r). 



isotherm are shown in Fig. [21 at densities of p* = 1.05 x 
10"'' and 1.00 x lO""'. The key point is that at these 
low densities, p{r) appears to be a superposition of two 
parts: a short-range associated-ion contribution, which 
dies off by about r = 2-3cr; and a peaked contribution 
corresponding to free ions. The short-range part decays 
within a distance much shorter than the Bjerrum length 
commonly used as a distance-based clustering criterion, 
which at this temperature is 20a. Assuming no correla- 
tions between free ions and any other ions in the system, 
the peaks in p{r) should occur at rp ~ 67a and 15a for 
p* = 1.05 X 10~^ and 1.00 x 10~^, respectively; by com- 
parison with the simulation results, these predictions are 
very reliable. The free-ion peaks are very broad, showing 
that a distance-based criterion rc ~ Ab is not physically 
justifiedji^ As the density of ions is increased, the peak 
both shifts to lower values of r and decreases in height. 
Equation ([TU|) provides excellent fits to the simulation 
results for p{r) (for r > rg), and yields values for the 
degree of association a as shown in Fig. [TJd). There 
is very good agreement with the simulation results us- 
ing Tc = 2(7, in correspondence with the comments made 
above regarding the decomposition of p{r) in to a short- 
ranged associated-ion contribution and a broad free-ion 
contribution. Accordingly, there is excellent agreement 
between the theoretical predictions for a and the results 
of the analysis of p{r). 

Figure [T] shows that the theoretical predictions are not 
very sensitive to the precise values of rc and hence K. 
This is explored further in Fig. |31 which shows K as 
a function of temperature evaluated using Eq. ^ with 
fc = Ab and Ab/2, and from the asymptotic expres- 
sion in Eq. ([7]). The first two expressions give essen- 
tially identical numerical results over the temperature 
range 0.04 < T* < 0.10, the region of current interest. 
The asymptotic expression, Eq. ([T]), is accurate only at 
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FIG. 2: Nearest-neighbour distribution function p(r) at T* — 
0.05. The simulation results are from simulations with p* = 
1.05 X lO"*^ (open symbols) and p* = 1.00 x 10~* (filled sym- 
bols); the thick curves are fits using the dissociated-ion result 
in Eq. (|10|) over the range r > ro, where ro is the maximum 
in p{r). 



the lower end of the temperature range. Note that K is 
plotted on a logarithmic scale; the deviation between the 
'Bjerrum length' and asymptotic results at T* = 0.07 is 
about 20%. 

When considering the effective interactions between 
charged colloids, it is of primary importance to know 
the degree of association of counterions and added salt 
within the suspending phasei^iii"— Of course, detailed 
calculations are easily performed using the prescriptions 
outlined herein. For a qualitative assessment, however, it 
is useful to divide the phase diagram into regions where 
the ions are mostly dissociated (a < ^) and where they 
are mostly associated (a > ^). Within the theory out- 
lined in Section ini the dividing line is defined by Eq. 
The phase diagram in the p*-T* plane is shown in Fi g. \M 
the vapour-liquid coexistence data are taken from Ref. l22l . 
The ct = ^ line is shown for the three expressions for K, 
with rc = Ab and Ab/2, and from Eq. ([7]). The devia- 
tions between these expressions become more pronounced 
as temperature increases, due to the increasingly signifi- 
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FIG. 3: Cation-anion configurational integral, K, as a func- 
tion of temperature T*: Eq. ^ with rc — Ab (solid line); 
Eq. (O with rc = Ab/2 (dashed line); Eq.Q (dot-dashed 
line) . 
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Salt concentration c (M) 
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FIG. 4: Phase diagram of the RPM showing the vapour-hquid 
coexistence region (lower right) from simulations,— and the 
boundary between dissociated (q < |) and associated (a > |) 
regimes as predicted from Eq. with three different evalu- 
ations of K: Eq. (|3} with Tc = Ab (black solid line) and 
= Ab/2 (red dotted line); Eq. (0 (green short-dashed line). 
Also shown is the locus of maxima in the const ant- volume heat 
capacity Cv from a simple two-particle theory^ (blue long- 
dashed line) and the boundary between ideal and strongly cor- 
related regimes predicted by DH theory from Eq. (|12|l (black 
dot-dashed line). Real units are shown for a monovalent salt 
with ionic diameter cr = 4 A at a temperature T = 298.15 K. 



cant large-r contributions to the integral in Eq. 

It should be noted that significant improvements on the 
non-interacting particle theory are possible, and indeed 
have been developed in detail. Already in 1926, Bjer- 
rum took account of the reduction in free-ion concentra- 
tion in order to compute the mean activity coefficients 
in electrolyte solutions Fisher and Levin have ex- 
plored the consequences of this level of approximation on 
the thermodynamics and phase behaviour of the RPM 
in their Debye-Hiickel-Bjerrum theory, and found that 
it leads to a vapour-liquid coexistence curve of the in- 
correct shapei ^^i^" Including ion-ion pair (ion-dipole) in- 
teractions restores the correct shape of the coexistence 
curve, and yields quite accurate values for the critical pa- 
rameters. These extensions are not applied here to the 
problem of ion pairing: as Figs. [Ijc) and (d) show, the 
agreement between theory and simulation at the low con- 
centrations and low temperatures of interest is excellent; 
the simplest ion-pairing theory is clearly adequate for the 
present purposes. 

The DH expression for the osmotic pressure of the elec- 
trolyte, in RPM reduced units, veadsS^'^^ 
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where x = = \/ 4:TTp* /T* and is the Debye 

screening length. This expression is exact to leading or- 
der in X as p* — ^ 0. At low density and high temperature, 
the ions are largely dissociated and the thermodynamics 
is essentially ideal. Significant deviations from ideality 
are expected when x/T* ^ 1, and a dividing line between 
frec-ion and strongly correlated regimes can therefore be 
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This line is included in Fig. and shows that there is 
a significant portion of the phase diagram in which the 
DH theory would suggest a low degree of ion association 
(because x/T* < 1), but the simulations and the Bjerrum 
theory show that a > ^. This has a serious consequence 
for linearised Poisson-Boltzmann theories of electrolyte 
solutions and related systems, which assume weak ion- 
ion correlations; in that part of the phase diagram lying 
between the lines defined by Eqs. ([6|) and (|T2|) . ions are 
associated and hence strongly correlated despite the fact 
that x/T* <1. 

There is one more feature of the a = ^ locus to be dis- 
cussed, and that is its monotonic variation with density. 
At high enough density, the distinction between two free 
ions and one ion pair becomes blurred. At the simplest 
level, this volume effect can be captured by a two-particle 
theory, in which all ions arc resolved into cation-anion 
pairs, and each ion pair has an internal configurational 
integral given by 



<?2 = 47r ; s 

la/2 



ds 



(13) 



where s is the distance from the ion-pair center of 
mass to one of the constituent ions, and the upper 
limit Sc = (3/27rp)^/'^ fixes the volume per pair to be 
2/ p. Although the degree of association is not defined 
within this theory, one can delineate the boundary be- 
tween dissociated and associated regimes with the lo- 
cus of maxima in the constant-volume heat capacity 
Cv = kBP^{d'^lnq2/dp'^)v, where /? = l/kuTM This 
line is shown in Fig. 2] and suggests that the domain of 
associated ions is bounded from above (T* ~ 0.1). This is 
in correspondence with the types of phase diagrams pro- 
posed for a wide range of ionic fluidsJ^ It was shown in 
earlier work that the two-particle theory provides an ex- 
cellent account of simulation measurements for the max- 
ima in Cv^ 



B. Analysis of experimental data 

To aid comparisons with experimental systems, axes 
in Fig. |4] are also shown with real units for an elec- 
trolyte with ionic diameter tr = 4 A at a temperature 
T = 298.15 K. The concentration Cg is given in moles 
of salt (not ions) per litre, and with the physical tem- 
perature held constant, T* becomes proportional to the 
solvent dielectric constant e. 

To effect a direct comparison between experiment and 
theory (and therefore to link experiment with simula- 
tion), attention is turned to the degree of association, 
a. This can be extracted from experimental measure- 
ments of the molar conductivity. Ignoring the formation 
of ion triples and higher charged clusters, a = 1 — A/Aq 
where A is the molar conductivity and Aq is its limiting 
value at infinite dilution. The most common source of 
experimental uncertainty is in the determination of Aq, 
since a suitable model has to be used to extrapolate A to 
infinite dilution. For the present purposes, conductivity 
data for tetraalkylammonium salts in various solvents are 
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FIG. 5; Degree of association a against Kp = 2KmCs/c^ for 
some tetraalkylammonium salts in various solvents. Refer- 
ences for the experimental data are given in square brack- 
ets. The solid line is Eq. ([Sj. (Bu =butyl, Pr =propyl, 
Pi =picrate.) 



analysed using the quoted estimates for Ag. In particu- 
lar, Leunissen and co-workers have obtained the degree of 
association of tetrabutylammonium halides in bromocy- 
clohexane and decalin-bromocyclohexane mixturesj ^"'^^ 
Conductivity data for tctrapropylammonium picrate in 
chlorobenzene^ and tetrabutylammonium iodide in car- 
bon tetrachloridc-nitrobenzenc mixtures^ arc also anal- 
ysed. The system parameters are summarised in Tabic H] 
the experimental data span wide ranges of concentration 
and solvent dielectric constant. 

The experimental results for a are fitted using Eq. ([5]), 
with the association constant as a fitting parameter. In 
fact, for the purposes of analysis, Kp is replaced by 
2KmCs/c^ , where Km is the dimensionless ion-pairing 
equilibrium constant (defined on the molar scale), 2cs is 
the ion concentration in M, and = 1 M. The fit- 
ted values of Km are reported in Table HI The experi- 
mental values for a are shown in Fig. [5l plotted against 
Kp = 2KmCs/c^ ■ There is an impressive collapse of the 
experimental data on to the theoretical universal curve 
given by Eq. 

It is possible to convert the experimental parameters 
in to RPM units. Since g, e, T, and hence Ab are all 
known, an effective ion diameter a can be obtained by 
equating the fitted values of Km with the expression for 
K in Eq. ([3]), the numerical conversion between the two 
being 



K 



K 



M 



lOOOiVAC^ 



(14) 



where Na is Avogadro's number. The integrand in 
Eq. (121) is exp (Ae/r), and Vc = Ab/2. (Other sen- 
sible choices for Tc give essentially the same numerical 
results.) The ion diameter a is obtained by numerical so- 
lution of Eq. This procedure ignores chemical detail 
such as the precise nature of the short-range repulsive 
interactions and the presence of van der Waals interac- 
tions, and so the fitted values of cr might reflect some 
non-coulombic effects. Nonetheless, under the conditions 
considered here, the coulombic interaction between ions 
at contact is dominant (being in the range of 10-20 fceT) 



and so the effective RPM parameters are expected to be 
meaningful. 

Values of Ab, cr, T* = ct/Ab, and the range of RPM 
ion densities p* are all summarised in Table HI The data 
of Leunissen and co-workers correspond to RPM temper- 
atures in the region of T* = 0.05, and RPM ion densities 
as low as p* ~ 10~*; Lindback and Beronius' datai^ cor- 
respond to a similar temperature, but do not extend to 
as low concentration. The data of Roy et al^ do not 
extend to very low concentration, but they do span the 
temperature range 0.06 < T* < 0.11. Note that the ex- 
perimental data correspond to the regime where the mean 
ion-ion separation {p~^^^) is much greater than the Bjer- 
rum length, i.e., in RPM units y^A^Fp* <C T* . This corre- 
sponds to the weakly interacting regime of the DH theory 
summarised in Section [IV A[ and so the theory would not 
have predicted the extensive pairing apparent in experi- 
ments. The central conclusion is that the Bjerrum theory 
(and by association, the simulations) successfully treats 
the ion-pairing equilibrium of nanomolar salt solutions in 
low-polarity solvents, such as those employed in recent 
work on chargcd-coUoid suspensions i^^— 

The effects of ion pairing on the effective screening 
length will now be considered. The appropriate effec- 
tive screening parameter in the presence of association is 
clearly 



A/47r(l - a)p\B = kdV^ 



(15) 



where kd = -s/^ttpAb is the inverse of the Dcbye length. 
At very low salt concentrations where Kp <C 1, Eq. ([S]) 
shows that a 



^Kp, and hence 
-^l-\Kp. 



(16) 



The screening length is then essentially equal to that as- 
suming complete dissociation, i.e., the Debye length. At 
higher concentrations where Kp ^ 1, Eq. ([S|) predicts 
a ~ 1 — ^/2/Kp and hence 



K 



2 
lip 



1/4 



(17) 



In this regime, k ^ kd and hence the screening length is 
much longer than would be expected on the basis of com- 
plete dissociation. Moreover, there is an unusual scaling 
behaviour: because kd has a p^^^ dependence, the effec- 
tive scaling parameter scales like k ~ p^^'^ ■ p~ 
Hence, the effective screening length scales like 



P 



1/4 



-1/4 



while the Debye length scales like Kj-,^ 



-1/2 



A glance at Fig. [5] shows that the regime Kp ^ 1 is in 
fact experimentally accessible. In particular, the exper- 
imental measurements by Leunissen and co-workersiSiS^ 
extend to a very high degree of association, and so the re- 
sulting interactions between charged colloids suspended 
in these solutions will not be screened as effectively as 
might be expected. As an example. Fig. [5] shows the De- 

T 1/2 

bye length Kj-, cx Cs , and the effective screening length 
from Eq. ((TS]), for a solution with cr = 4 A, e = 7, 
and T = 298.15 K; this set of parameters corresponds to 
a reduced RPM temperature of T* = 0.05, and is rep- 
resentative of experiments in low-polarity solvents. Over 
the range 10~^ M < Cs < 10~^ M, the screening length 
and the Debye length coincide, ranging from several mi- 
crometres down to several hundred nanometres; these 
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TABLE I: Physical parameters for tetraalkylammonium salts in various solvents: e is the solvent dielectric constant, Ab is the 
Bjerrum length for T = 298.15 K, Cs is the salt concentration, A'm is the dimensionless ion-pairing association constant (defined 
on the molar scale), and a is the effective hard-sphere diameter of the ion. T* = (t/Ab and p* are the effective RPM temperature 
and density, respectively. (Bu =butyl, Pr ^propyl, Pi =picrate.) 



System 






Ref. 


e 


Ab (A) 


Cs i/iM) 


Km 


a (A) 


T* 


lOV* 


BU4N+CI 


- / CsHiiBr 




10 


7.92 


70.8 


0.28-244 


3.80 X 10^ 


3.27 


0.0462 


0.0118-10.3 


Bu4N+Br 


- / CeHiiBr 




35 


7.92 


70.8 


0.56-245 


9.17 X 10® 


3.57 


0.0504 0.0307-13.4 


Bu4N+Br 


- / 27.3 wt% decalin-CeHiiBr 


35 


5.62 


99.7 


0.47-157 


6.31 X lO'' 


5.53 


0.0555 


0.0957-32.0 


Pr4N+Pi- 


/ PhCl 




44 


5.612 


99.9 


5.97-1539 


8.80 X 10*^ 


5.41 


0.0542 


1.14-294 


BU4N+I- 


/ 80 wt% CCI4 


-PhN02 


45 


10.22 


54.8 


170-1330 


6.44 X 10^ 


3.16 


0.0576 


6.46-50.5 


BU4N+I- 


/ 60 wt% CCI4 


-PhN02 


45 


17.45 


32.1 


170-1320 


1.23 X 10^ 


3.17 


0.0987 


6.52-50.6 


BU4N+I- 


/ 40 wt% CCI4 


-PhN02 


45 


23.90 


23.5 


440-4070 


2.91 X 10^ 


2.60 


0.111 


9.31-86.2 


BU4N+I- 


/ 20 wt% CCI4 


-PhN02 


45 


29.66 


18.9 


830-7380 


1.57 X 10^ 


2.08 


0.110 


9.00-80.0 



values correspond well with experimentally determined 
values.—!^ At higher salt concentrations Cg > 10^^ M, 
exceeds due to the formation of ion pairs. In 
addition, decays less fast with increasing concentra- 
tion, switching over to the Cs dependence advertised 
in Eq. (|17|) . Finally, it is noted that at very high con- 
centrations (outwith the relevant range studied here) the 
formation of ion pairs can lead to increasing with 
increasing Cg; the ion-pair contribution to the effective 
dielectric constant of the solution increases with increas- 
ing concentration, which ultimately leads to reductions 
in Ab and kJ^ 

Some interesting transient behaviour has been observed 
in experiments on chargcd-colloid suspensionsi^ Specifi- 
cally, the colloids show the effects of anomalously long- 
ranged repulsions which are incompatible with the ap- 
parent salt concentration Cs ~ 10~^ M. This transient 
behaviour can occur on the timescale of a few days. One 
contributing factor might be the time taken for ions to as- 
sociate and/or dissociate after preparation. For instance, 
if ions were initially associated, then the screening length 
and colloidal repulsions would decrease on the approach 
to equilibrium. Assuming that the association of cations 
with anions is a second-order, diffusion-controlled 'reac- 
tion', then the corresponding macroscopic rate constant 
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FIG. 6: The Debye length k.^^ , and the effective screening 
length from Eq. (|15p . for a solution with ct = 4 A, e = 7, 
and T = 298.15 K; this set of parameters corresponds to a 
reduced RPM temperature of T* = 0.05. 



can be estimated by fca = 8RT/3r]^ where 77 is the vis- 
cosity of the solvent. For a solvent with rj = 10~^ Pa s 
at room temperature, this yields fca = 7 x 10^ M^^ s^^. 
The corresponding first-order dissociation of ion pairs will 
have a rate constant fed = ^aC^/i^M- Table U shows 
that Km can be as high as 10^, and so the charac- 
teristic time for ion dissociation will be no more than 
1/fcd ~ 0.001 s. This is only a very rough estimate and 
solvation-shell structure and the slow escape from the 
long-range coulombic attraction between ions may well 
increase this timescale considerably, but it seems unlikely 
that ion dissociation is the dominant cause of the tran- 
sient behaviour. A more prosaic explanation for the ob- 
served transients in Ref. y may be that ions are initially 
sequestered but subsequently released by impurities such 
as water or by the container walls, leading to a slow re- 
duction in the screening length. 



V. CONCLUSIONS 

In this paper, the association of ions in low-polarity 
solvents was studied within the context of the restricted 
primitive model at low concentration and low tempera- 
ture. In reduced units, concentrations as low as 10~^° 
and temperatures as low as 0.04 were simulated using an 
efficient Monte Carlo algorithm^^i and comparisons were 
made with Bjcrrum's quasi-chemical ion-pairing theoryji^ 
These conditions correspond to nanomolar salt solutions 
in low dielectric constant solvents under ambient condi- 
tions; such media are of relevance to recently synthesised 
charged-coUoid suspensions 1^"— The degree of ion asso- 
ciation strongly affects the effective interactions between 
colloids r^"— and so one of the aims of this work was to 
map out the different regimes of association on the phase 
diagram. 

In the simulations, distance-based criteria and the 
nearest-neighbour cation-anion distribution function p{r) 
were used to determine the degree of association. It was 
shown that conventional distance-based criteria based on 
the Bjerrum lengtbi^ are inferior to a short-range cut- 
off J^i^ This conclusion was backed up by analysis oip{r), 
which shows that associated ions are in close proximity, 
while free ions show a very broad distribution of distances 
to the nearest ion of opposite charge. 

To some extent, this feature is accounted for in the- 
oretical treatments based on an ion-pair configurational 
integral, since it is the Boltzmann factor at short-range 
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which makes the most significant contribution, at least 
at low temperatures where association is prevalent. The 
precise value of the upper limit in the integral is not im- 
portant, as already noted long agoi^ The agreement be- 
tween Bjcrrum theory and simulation under the physi- 
cal conditions studied here is excellent. More sophisti- 
cated treatments, such as those developed by Fisher and 
Levin^iSS are scarcely required; of course, near to the 
coexistence region, ion-ion and ion-dipole interactions be- 
come of paramount importance. Additionally, the anal- 
ysis of experimental dielectric and conductivity data at 
low temperature and moderate concentrations requires 
a detailed account of ion triples and higher (charged) 
clusters?^ such phenomena as conductivity minima in 
the vicinity of the coexistence region can be handled us- 
ing suitable extensions of the DH theory that include the 
effects of association42i^ 

A comparison of the Bjerrum theory and the Debye- 
Hiickel theory shows that there is a significant region 
of the phase diagram where the former predicts strong 
ion association (in agreement with simulations) but the 
latter indicates only weak ion-ion correlations. This im- 
plies that linearised Poisson-Boltzmann theories, which 
assume weak ion-ion correlations, should only be ap- 
plied under conditions where ion association is actually 
known to be insignificant. Those conditions are identi- 
fied in this work. Of course, these restrictions should also 
be observed in the application of the Derjaguin-Landau- 
Verwey-Overbeek theory because it incorporates Debye 
screening of colloidal interactions under the assumption 
of weak ion-ion correlations, at least for the case of spher- 
ical colloidal particles. 

The Bjerrum theory was compared with experimental 



data for tetraalkylammonium salts in low-polarity sol- 
vents, and the agreement was found to be excellent. Fit- 
ting the association constants allowed a mapping between 
real systems and the restricted primitive model. Simula- 
tion, theory, and experiment have therefore been com- 
pared at very low concentrations and low temperatures. 
One of the primary motivations for this study was to un- 
derstand the association of ions in low-polarity solvents, 
with a view to getting a clear picture of the nature of 
screening between colloidal particles suspended in such 
solutions. The effects of ion association on the appro- 
priate screening lengths have been quantified, and the 
results are in complete accord with experiment. It is 
hoped that the results presented herein can be used in 
subsequent treatments of charged-coUoid interactions in 
low-polarity media. 
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Figure captions 

Fig. 1 The degree of association a along isotherms. The curves are the theoretical predictions of Eq. ([S]) using K 
computed with a cut-off Tc = Ab/2 (soUd lines), and using the asymptotic result in Eq. (O for K (dashed 
lines). The points are the simulation results computed using various criteria: (a) distance criterion with 
Tc — Ab; (b) distance criterion with Tc = Ab/2; (c) distance criterion with Tc = 2a; (d) fitting Eq. (|T0)) to 
p{r). 

Fig. 2 Nearest-neighbour distribution function p(r) at T* = 0.05. The simulation results are from simulations with 
p* ~ 1.05 X 10~^ (open symbols) and p* = 1.00 x 10^* (filled symbols); the thick curves are fits using the 
dissociated- ion result in Eq. (|10p over the range r > ro, where rg is the maximum in p{r). 

Fig. 3 Cation-anion configurational integral, K, as a function of temperature T*: Eq. ^ with rc = Ab (solid 
line); Eq. ([3]) with rc = Ab/2 (dashed Hne) (this is almost indistinguishable from the former curve, in this 
temperature range); Eq.([7]) (dot-dashed line). 

Fig. 4 Phase diagram of the RPM showing the vapour-liquid coexistence region (lower right) from simulations ^2 
and the boundary between dissociated (a < ^) and associated (a > ^) regimes as predicted from Eq. © 
with three different evaluations of K: Eq. ([3]) with = Ab (black solid line) and rc = Ab/2 (red dotted 
line); Eq. ([7]) (green short-dashed line). Also shown is the locus of maxima in the constant- volume heat 
capacity Cy from a simple two-particle theorj*^ (blue long-dashed line) and the boundary between ideal 
and strongly correlated regimes predicted by DH theory from Eq. (fT^ (black dot-dashed line). Real units 
are shown for a monovalent salt with ionic diameter cr = 4 A at a temperature T = 298.15 K. 

Fig. 5 Degree of association a against Kp = 2KmCs/c'^ for some tetraalkylammonium salts in various solvents. 

References for the experimental data are given in square brackets. The solid line is Eq. ([5]). (Bu =butyl, 
Pr ^propyl. Pi =:picrate.) 

Fig. 6 The Debye length n^^, and the effective screening length k^^ from Eq. ([T5|). for a solution with tr = 4 A, 
e = 7, and T = 298.15 K; this set of parameters corresponds to a reduced RPM temperature of T* = 0.05. 



